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I.  INTRODUCTION  AND  SUMMARY 
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1.1  BACKGROUND 

The  objective  of  the  Systems,  Science  and  Software  (S3) 
research  program  is  to  examine  the  parameters  that  affect  the 
seismic  signals  from  underground  explosions  and  earthquakes. 
Attention  is  primarily  directed  to  those  features  of  the 
seismic  waveforms  that  discriminate  between  the  two  classes 
of  events  and  that  reliably  indicate  the  explosion  yield. 
Current  research  includes  empirical  studies  of  the  available 
data,  time  signal  analysis,  and  the  development  and  application 
of  theoretical  and  numerical  methods  for  modeling  earthquakes 
and  explosions.  Emphasis  is  on  the  last  of  these.  In  partic- 
ular, we  are  applying  techniques  for  theoretically  simulating 
the  far-field  signatures  of  simple  and  complex  seismic  sources. 

This  report  summarizes  the  work  done  during  the  seventh 
three-month  period  of  the  current  contract. 


1.2  SUMMARY  OP  RESEARCH  DURING  THIS  QUARTER 

Our  work  during  this  quarter  has  included  research  in 
a number  of  areas.  Research  projects  currently  underway  or 
completed  during  the  quarter  are  briefly  summarized  below. 

Source  Studies 

A.  Earthquake  Modeling  on  the  ILLIAC  IV  Computer 

A three-dimensional  finite  difference  program  for 
modeling  earthquake  faulting  has  been  made  operational  on 
the  ILLIAC  IV  computer.  The  current  version  is  capable  of 
handling  a bilateral  fault  in  a homogeneous  medium.  Our 
initial  set  of  calculations  is  designed  to  understand  the 
importance  of  plastic  material  behavior  in  the  fault  zone 
and  to  determine  the  scaling  of  the  radiation  field  with 
the  fault  parameters. 
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Two  calculations  have  been  successfully  completed.  Both 
are  for  a 3 km  x 3 km  bilateral  fault  in  a wholespace.  In  one 
case  the  material  behavior  is  linear  elastic.  In  the  other 
case  plastic  yielding  is  allowed  in  the  vicinity  of  the  fault 
plane. 

Detailed  analysis  of  the  linear  elastic  calculation  has 
been  completed  and  excellent  results  were  obtained.  The  far- 
field  radiation  was  computed  with  two  methods.  One  method  is 
to  use  dislocation  theory  with  the  computed  slip  time  histories 
on  the  fault.  The  second  method  is  to  expand  the  outgoing  wave 
field  in  spherical  harmonics  (multiple  coefficients) . Only  the 
latter  is  formally  correct  when  nonlinear  material  behavior  is 
allowed  near  the  fault. 

The  second  (nonlinear)  earthquake  simulation  is  now 
being  analyzed  the  same  way  as  the  first  and  results  from  the 
two  are  being  compared.  A detailed  report  describing  this  work 
is  being  written. 

B.  Representation  Theorem  for  Analytic  Continuation  of 

Finite  Difference  Source  Calculations 

This  work  is  summarized  in  Section  1.4  and  described 
in  detail  in  Section  III. 

c.  Source  Calculations 

This  work  is  summarized  in  Section  1.5  and  described 
in  detail  in  Section  IV. 


Data  Analysis 

A.  Discrimination  Experiment 

This  work  is  summarized  in  Section  1.3  and  described 


in  detail  in  Section  II. 


The  Variable  Frequency  Magnitude  (VFM)  discriminant 
being  applied  in  the  discrimination  experiment  (see  Section  II) 
is  based  on  the  use  of  narrow  band  filter  magnitudes,  called 
m^Cf).  These  magnitudes  reflect  the  amplitude  of  the  arriving 
energy  at  a particular  period  and  a particular  arrival  time. 

In  essence,  they  represent  the  spectral  energy  of  a specific 
phase  arrival.  We  have  long  recognized  that  measurements  of 
this  kind  might  provide  a more  convenient,  stable  and  reliable 
signal  amplitude  indicator  than  the  conventional  m^. 

A set  of  HNME  recordings  of  eleven  Pahute  Mesa  recordings 
was  provided  in  digital  form  and  is  being  examined  to  test  the 
utility  of  the  m^(f) . A similar  value,  Mg(f),  is  defined  for 
the  surface  waves.  Single  frequency  values  are  no  improvement 
on  the  conventional  time  domain  magnitudes.  However,  m^Cf)  and 
Mg  (f ) values  defined  by  averaging  the  m^f)  and  Mg(f)  over  a 
frequency  band  are  very  attractive  in  terms  of  such  qualities 
as  the  scatter  when  plotted  versus  yield.  We  are  continuing  to 
evaluate  this  measurement  and  will  report  the  results  to  VSC  in 
the  near  future. 

Surface  Waves 

A.  Surface  Wave  Amplitudes  of  NTS  Explosions  Recorded  at 

ALQ  and  TUC 

A special  topical  report  describing  this  work  has  been 
submitted  to  VSC.  The  abstract  of  this  report  is  reprinted 
below.  A classified  report  with  detailed  tabulations  of  the 
important  results  has  also  been  prepared  and  will  be  submitted 
to  VSC  in  the  near  future. 

"Source  Amplitudes  of  NTS  Explosions  Inferred  from 
Rayleigh  Waves  at  Albuquerque  and  Tucson,"  by  T.  C.  Bache, 

W.  L.  Rodi  and  B.  F.  Mason. 


Abstract 


Comparing  observed  and  synthetic  seismograms,  source 
amplitudes  of  NTS  explosions  are  inferred  from  Rayleigh  wave 
recordings  from  the  WWSSN  stations  at  Albuquerque,  New  Mexico 
(ALQ)  and  Tucson,  Arizona  (TUC) . The  potential  influence  of 
source  complexities,  particularly  surface  spallation  and 
related  phenomena,  is  studied  in  detail. 

As  described  in  earlier  work  by  Bache,  Rodi  and  Harkrider, 
the  earth  models  used  in  computing  the  synthetic  seismograms 
were  inverted  from  observations  at  ALQ  and  TUC.  The  agreement 
of  observed  and  synthetic  seismograms  is  quite  good  and  is  sensi- 
tive to  important  features  of  the  source. 

The  events  studied  are  in  three  distinct  areas.  Yucca 
Flat,  Pahute  Mesa  and  PILEDRIVER  in  Climax  Stock  granite.  All 
events  were  below  the  water  table  and  the  yields  varied  from  40  to 
200  KT.  Within  each  group  the  mean  static  value  of  the  reduced 
displacement  potential  was  determined  at  a fixed  scaled 

yield,  assuming  a spherically  symmetric  point  source.  The 
source  is  then  modified  to  study  the  effect  of:  (1)  the  addi- 
tion of  a do uble- couple  component;  (2)  the  addition  of  a sur- 
face impulse  associated  with  impact  of  the  material  spalled 
from  the  surface;  (3)  loss  of  energy  from  the  free  surface 
reflected  waves  due  to  spallation.  Comparing  observed  and 
synthetic  seismograms , the  extent  of  the  effect  of  these 
secondary  phenomena  is  outlined. 

For  Yucca  Flat  and  Pahute  Mesa  events  the  inferred  source 
amplitudes  are  in  general  agreement  with  values  obtained  by 
other  methods.  The  spall  impulse  is  too  small  to  be  of  much 
importance.  More  likely  to  be  important,  especially  for  the 
Pahute  Mesa  events,  is  the  loss  of  energy  from  the  upgoing  waves. 
For  PILEDRIVER  the  double-couple  and  attenuation  of  upgoing 
waves  dominate  the  source  determination.  Correcting  for  these 
phenomena,  the  explosion  source  level  (4^)  is  considerably 
smaller  than  has  usually  been  supposed. 


Body  Wave  Studies 

In  the  section  on  "Data  Analysis"  we  mentioned  the  set 
of  eleven  HNME  recordings  of  Pahute  Mesa  explosions  that  are 
being  used  to  test  the  development  of  newly  defined  magnitude 
measures.  These  data  are  also  being  analyzed  with  synthetic 
seismogram  methods  to  isolate  the  different  effects  contri- 
buting to  mb-log  yield  and  Ms~log  yield  relations.  An  impor- 
tant feature  of  the  body  wave  recordings  is  a depth-dependent 
phase  that  arrives  after  pP.  This  may  be  due  to  spall  slap- 
down,  tectonic  release  or  some  other  cause. 

To  facilitate  analysis  of  these  seismograms,  our  body 
wave  synthesis  code  (GENSRC)  has  been  improved.  The  current 
version  will  compute  the  far-field  body  waves  for  a source 
consisting  of  an  explosion  (specified  by  a reduced  displace- 
ment potential)  plus  a tectonic  release  component  (specified 
by  a double-couple)  plus  a spall  slapdown  component  (specified 
by  a stress-time  history  applied  to  the  free  surface)  in  a 
layered  half space. 

1.3  SUMMARY  OF  SECTION  II:  "DISCRIMINATION  EXPERIMENT" 

During  this  reporting  period  we  received  short-period 
seismograms  in  digital  format  for  nine  additional  Eurasian 
events,  bringing  the  total  number  of  events  in  the  data  set 
to  28.  The  event  identification  study,  based  on  variable 
frequency  magnitude  estimates,  was  expanded  to  include  eight 
of  the  stations  that  have  contributed  the  bulk  of  the  data. 

In  addition,  noise  corrections  were  applied  to  the  magnitude 
estimates  at  one  of  the  stations  (Bluff,  Alaska)  and  this 
resulted  in  a significant  enhancement  in  the  separation  of 
the  earthquake  and  explosion-like  populations.  Tentative 
identifications  for  the  28  events  based  on  the  combined  re- 
sults from  the  stations  analyzed  are  as  follows:  earthquakes  - 
events  36,  38,  41,  47,  48,  49,  50,  55  through  62  and  64; 


explosion-like  1,  14,  16  through  22,  33,  53;  event  63  - mixed 
identification  (possibly  deep) . The  event  locations  are 
described  in  Figure  2.1  and  Table  2.1  of  Section  II.  Future 
analysis  will  include  noise  corrections  at  all  the  stations 
and  multi-station  discrimination. 

1.4  SUMMARY  OF  SECTION  III:  "SURFACE  WAVE  SOLUTION  FOR  THE 

ELASTIC  EQUIVALENT  OF  A COMPLEX  AXI SYMMETRIC  SOURCE" 

An  important  problem  in  computing  theoretical  seismograms 
for  complex  deterministic  source  models  is  the  linking  of  finite 
difference  numerical  source  models  with  analytical  methods  for 
propagating  elastic  waves  in  realistic  earth  models.  In  Section 
? Ill  we  outline  the  mathematical  development  of  a method  to 

analytically  continue  the  stress  waves  from  an  axisymmetric 
source  into  a layered  earth  model  for  computation  of  far-field 
surface  waves.  The  method  presented  allows  for  computation  of 
Rayleigh  wave  motion  given  the  Fourier  transformed  displacements 
and  stresses  monitored  on  the  edge  of  a cylinder  surrounding  the 
source  region.  Though  the  resulting  expressions  appear  compli- 
cated, they  have  an  asymptotic  form  which  should  make  interpre- 
tation straightforward.  The  surface  wave  methods  presented  here 
together  with  the  body  wave  methods  presented  in  an  earlier 
report  will  allow  efficient  computation  of  theoretical  seismo- 
grams at  nearly  all  distances  of  interest. 

1.5  SUMMARY  OF  SECTION  IV:  "SOURCE  CALCULATIONS" 

As  a prelude  to  our  two-dimensional  source  calculations, 
we  have  initiated  a review  of  our  one-dimensional  spherically 
symmetric  source  calculations  in  order  to  update  and  improve 
our  nonlinear  constitutive  models  while  ensuring  that  they  are 
compatible  with  our  two-dimensional  codes.  The  three  most 
important  modeling  improvements  are:  (1)  the  implementation 
of  good  equations  of  state  for  the  cavity  gases  for  salt, 
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granite  and  tuff;  (2)  improvements  which  make  the  overburden 
pressure  treatment  consistent  with  our  equation  of  state,  and 
(3)  the  development  of  a better  effective  stress  law  to  model 
the  influence  of  water  saturation  in  rocks.  In  Section  IV, 
we  present  discussions  of  our  latest  one-dimensional  calcula- 
tions for  PILEDRIVER  (granite) , SALMON  (salt)  and  Rainier  Mesa 
saturated  tuff,  and  discuss  the  modeling  improvements  made  for 
each. 

We  are  unable  to  match  the  SALMON  reduced  displacement 
potential  (RDP)  with  the  new  models  using  laboratory  data  for 
the  failure  envelope.  The  calculations  give  the  measured 
cavity  radius,  but  too  small  an  RDP.  We  will  now  attempt  to 
calculate  this  event  in  two-dimensions  to  study  the  influence 
of  in  situ  stress  in  the  RDP.  If  in  situ  stress  proves  not  to 
be  the  explanation  of  our  low  calculated  RDP,  we  plan  to  include 
viscoelastic  effects  in  the  constitutive  model  for  salt.  Our 
calculations  for  PILEDRIVER  and  for  saturated  tuff  were  in 
good  agreement  with  both  near  and  far  field  data. 
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II.  DISCRIMINATION  EXPERIMENT 


2 . 1 INTRODUCTION 

Our  objective  in  the  discrimination  experiment  is  to 
analyze  seismic  waveforms  from  a large  population  of  events 
in  order  to  identify  these  events  as  either  earthquakes  or 
explosions.  The  waveforms  that  we  are  concentrating  on  are 
short-period  P waves  recorded  by  a global  network  of  seismo- 
graph stations.  In  this  section  of  the  report  we  will  sum- 
marize the  work  that  has  been  performed  to  date. 


2.2  DESCRIPTION  OF  DATA 

As  of  the  end  of  this  reporting  period  we  had  received 
short-period  digitally  recorded  seismograms  for  28  Eurasian 
events  recorded  at  one  or  more  of  a network  of  20  globally 
distributed  seismograph  stations.  The  event  locations  are 
indicated  in  Figure  2.1  by  the  solid  circles  and  the  event 
dates,  origin  times  and  epicentral  coordinates,  as  supplied 
by  Teledyne  Geotech,  are  listed  in  Table  2.1.  With  the 
increase  in  the  number  of  events  in  the  data  base,  several 
concentrations  of  activity  have  become  apparent.  The  most 
active  region  is  the  Kuril  and  Kamchatka  Islands  which 
accounts  for  nearly  one-half  of  the  available  events.  The 
locations  of  the  eight  Eurasian  seismograph  stations  providing 
data  for  this  experiment  are  indicated  by  the  solid  triangles 
in  Figure  2.1. 

The  status  of  the  data  base  as  of  June  30,  1978  is 
described  in  Table  2.2  on  a station-by-station  basis  for  each 
of  the  28  events.  As  can  be  seen,  the  number  of  station 
seismograms  is  highly  variable  from  event  to  event,  ranging 
from  two  stations  for  event  36  to  twelve  stations  for  events 
20  and  53.  To  date  no  short-period  seismograms  have  been 


Figure  2.1.  Map  of  Eurasia  showing  locations  of  events  (solid  circles)  and  several  of 
the  stations  (solid  triangles)  included  in  the  discrimination  experiment. 


TABLE  2.1 


EVENT  LOCATION  INFORMATION 


Coordinates 


Event 

Number 

Date 

Yr  Mo  Dy 

Origin  Time 

Hr :Mn:Sec 

Latitude 

(°N) 

Longitude 

(°E) 

1 

77  07  26 

16:59:59.9 

69.4 

90.4 

14 

77  07  30 

01:56:59.9 

49.7 

78.2 

16 

77  08  10 

22:00:00. 7 

50.9 

111.0 

17 

77  08  17 

04:26:59. 8 

49.8 

78.2 

18 

77  08  20 

22:00:00.6 

64.1 

99.8 

19 

77  09  01 

03:00:00 

73.0 

54.0 

20 

77  09  05 

03:03:00 

50.0 

78.9 

21 

77  09  10 

16:00:00 

57.0 

106.0 

22 

77  09  30 

06:59:00 

48.0 

48.0 

33 

77  10  09 

11:00:00 

73.0 

55.0 

36* 

77  10  16 

20:03:35 

48.4 

152.9 

38 

77  10  16 

15:02:49 

36.9 

71.5 

41 

77  10  13 

20:38:42 

38.1 

72.8 

47 

77  10  16 

21:05:35 

49.7 

155.1 

48 

77  10  19 

05:02:00 

36.3 

71.3 

49 

77  10  19 

21:20:37 

49.5 

155.4 

50 

77  10  20 

08:18:04 

56.3 

164.1 

53 

77  10  29 

03:07:00 

49.0 

78.0 

55 

77  10  26 

05:38:52 

49.0 

155.8 

56 

77  10  26 

07:11:31.3 

46.4 

153.5 

57 

77  10  26 

13:14:30.9 

51.5 

153.4 

58 

77  10  27 

07:20:28.9 

53.5 

160.0 

59 

77  10  28 

21:15:11 

39.8 

71.9 

60 

77  10  29 

04:14:56.5 

47.0 

152.3 

61 

77  10  29 

06:26:42.5 

41.0 

63.7 

62 

77  10  29 

10:33:59.4 

47.3 

153.1 

63 

77  10  30 

21:38:15.6 

44.8 

145.0 

64 

77  10  31 

09:40:03.5 

55.8 

162.7 

* 

Alaska 

Data  Only 

STATUS  OF  DATA  BASE:  JUNE  30,  1978 


received  from  the  stations  BOCO  and  SNZO,  SRO  stations  in  Columbia 
and  New  Zealand,  respectively.  In  addition,  while  seismograms 
from  the  station  KSRS  are  available  for  27  of  the  28  events  in 
the  data  base,  the  extremely  high  level  of  background  noise 
prevailing  at  that  site  has  resulted  in  seismograms  that  are 
dominated  by  noise  for  the  majority  of  events.  In  the  fol- 
lowing subsection  of  this  report  we  will  describe  the  discrim- 
ination results  obtained  to  date  at  several  of  the  stations. 


2.3  DISCRIMINATION  RESULTS 

As  described  in  the  last  quarterly  report  written  under 
this  contract  (Rodi,  et  al . , 1978),  we  employ  a narrow  band 
filtering  procedure  to  compute  estimates  of  body  wave  magni- 
tudes, ( f ) , at  several  different  frequencies  within  the 

teleseismic  bandpass  (e.g.,  0.3  to  several  Hertz).  By  com- 
paring low  frequency  magnitude  estimates  with  high  frequency 
estimates  we  can  test  for  event  discrimination  using  short- 
period  P waves  recorded  at  the  different  stations. 

Figure  2.2a  shows  a seismogram  for  a short-period  P 
wave  from  event  47  recorded  at  the  station  in  Bluff,  Alaska 
(BFAK) . The  signal  is  preceded  by  approximately  30  seconds 
of  background  noise.  The  key  feature  of  the  signal  analysis 
procedure  we  employ  is  the  use  of  narrow  band  filters  to 
decompose  a time  series  consisting  of  signal  plus  noise,  such 
as  that  in  Figure  2.2a,  into  a set  of  quasi-harmonic  modulated 
signals.  The  modulation  or  envelope  function  is  calculated  by 
means  of  the  Hilbert  transform  with  the  maximum  of  the  envelope 
function  occurring  at  the  group  arrival  time  (tg)  of  energy  at 
the  center  frequency  (f  ) of  a particular  filter.  The  ampli- 

C 

tude  of  the  envelope  function  is  proportional  to  the  spectral 
amplitude  of  the  filtered  signal.  The  Gaussian-shaped  narrow 
band  filters  ensure  optimal  time  and  frequency  domain  resolu- 
tion within  the  constraints  imposed  by  the  sampling  theorem  or 
uncertainty  principle. 


Short-period  seismogram  for  the  P wave  recorded  at  Bluff,  Alaska  (BFAK) 
from  event  47. 


Figure  2.2b  shows  the  tg  versus  f representation  of 
the  time  series  in  Figure  2.2a.  The  envelope  peaks  from  20 
narrow  band  filters,  ranging  from  0.4  Hz  to  5.0  Hz,  are  scaled 
according  to  their  relative  amplitudes.  For  instance,  in  the 
case  of  the  filter  center  frequency  f = 2.75  Hz  the  scheme 
for  labeling  the  peaks  occurring  between  tg  = 402.7  and  406.9 
seconds  (after  event  origin  time)  is  that  * corresponds  to  the 
maximum  amplitude,  9 to  an  amplitude  between  90  and  < 100  per- 
cent of  * and  8 to  an  amplitude  between  80  and  < 90  percent  of 
*.  The  arrival  of  an  undispersed  P wave  signal  will  appear  as 
a horizontal  alignment  of  large  amplitude  peaks  in  this  plane. 
The  prominence  of  the  arrival  depends  on  the  ratio  of  the 
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signal-to-noise  spectral  amplitude  over  the  frequency  band  of 
interest.  In  Figure  2.2b  an  undispersed  arrival,  corresponding 
to  the  P wave  in  Figure  2.2a,  is  observed  at  a tg  of  approxi- 
mately 403  seconds  over  the  frequency  band  0.5  Hz  to  2.75  Hz. 
Note  that  the  uncertainty  in  the  t^  estimates,  Atg,  is  given 
by  Atg  Q /j"c>  where  Af/f  = Q 1,  Af  is  the  half-width  of 

the  narrow  band  filter  at  half  power  and  Q = 10.  Thus,  the 
uncertainty  in  the  group  arrival  times  Atg,  increases  with 
decreasing  filter  center  frequency. 

The  envelope  peaks  occurring  at  times  earlier  than 
— P — P 

tg  - aAtg,  where  tg  is  an  average  group  arrival  time  for  the 
P wave  computed  over  a frequency  band  corresponding  to  an 
acceptable  signal-to-noise  ratio  and  a is  a parameter  > 1,  give 
estimates  of  the  background  noise  and  can  be  used  to  "correct" 
the  signal  peaks  used  for  fi^Cf)  discrimination  tests.  While 
not  evident  in  Figure  2.2b,  later  arrivals  (tg  > t g + aAtg) 
could  be  identified  with  this  procedure.  In  Figure  2.2c  the 
sum  of  the  amplitudes  of  the  envelope  peaks  at  the  20  different 
frequencies  is  plotted  as  a function  of  time.  This  plot  indi- 
cates the  ratio  of  signal  amplitude  to  the  noise  or  later 
arrivals. 
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Figure  2.2b.  Group  arrival  time  (tg)  as  a function  of  filter  center  fre- 
quency (fc)  for  event  47  recorded  at  Bluff,  Alaska  (BFAK) . 
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Figure  2.2c.  Sum  of  the  envelope  amplitudes  plotted  as  a function  of  time 

after  origin  time  for  event  47  recorded  at  Bluff,  Alaska  (BFAK) 


Figures  2.3a,  2.3b  and  2.3c  are  a similar  sequence,  that 

is  P wave  seismogram,  t -f  and  sum  of  envelope  peaks,  for  an 

g c 

event  that  was  previously  identified  (Rodi,  et  al.,  1978)  as 
being  explosion- like.  The  sigr.al-to-noise  ratio  for  this 
event  as  recorded  at  BFAK  is  much  larger  than  that  for  event 
47  (e.g.,  compare  Figures  2.2a  and  2.3a).  As  a result,  in 
Figure  2.3b  there  is  a near  perfect  horizontal  alignment  of 
tg  estimates  over  the  entire  frequency  band,  0.4  Hz  to  5.0  Hz, 
with  a mean  tg  = 616  seconds.  In  Figure  2.3c  the  sum  of  the 
envelope  peaks  indicates  a signal-to-noise  ratio  for  this 
event  of  approximately  45. 

The  procedure  for  testing  for  discrimination  in  the 
m^Cf)  plane  follows  from  the  above  results.  Magnitude  esti- 
mates at  several  combinations  corresponding  to  peaks  with  tg's 
within  "tg  + aAt  are  selected  for  the  different  events  and 
compared  on  a station-by-station  basis.  Plots  of  m^ff)  planes 
for  eight  different  stations  are  shown  in  the  following  figures. 

In  each  of  these  figures  the  arrows  attached  to  the  closed 
circles  or  triangles  indicate  the  principal  direction  that  a 
noise  correction  would  move  the  points.  The  length  of  an  arrow 
is  proportional  to  this  correction. 

Figures  2.4a  through  2.4c  shows  rn^Cf)  results  for  18  events 
recorded  at  Bluff,  Alaska  (BFAK) . The  epicentral  distances 
for  these  events  ranges  from  25  to  65  degrees.  Three  different 
combinations  of  high  and  low  filter  center  frequencies  are 
plotted  in  order  to  give  an  indication  of  the  behavior  of  the 
event  populations.  In  the  previous  quarterly  report  events 
1,  14,  16,  17,  21,  33  and  53  were  identified  as  possible  explo- 
sions. Thus,  the  lines  drawn  in  these  figures  mark  the  bound- 
aries between  earthquake  and  explosion- like  events.  Comparing 
these  three  figures  we  see  very  little  movement  of  all  the 
events  except  63,  and  to  a less  extent  events  64  and  16,  with 
no  crossings  of  the  population  boundaries.  The  other  point  to 
note  about  the  events  in  all  three  of  these  figures  is  that  in 
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Short-period  seismogram  of  P wave  recorded  at  Bluff,  Alaska  (BFAK)  from 
event  53. 
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Figure  2.3b.  Group  arrival  time  (tg)  as  a function  of  filter  cente 
frequency  (f  ) for  event  53  recorded  at  Bluff,  Alaska 
(BFAK) . 
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Figure  2.3c.  Sum  of  the  envelope  amplitudes  plotted  as  a function  of  time 
after  origin  time  for  event  53  recorded  at  Bluff,  Alaska 
(BFAK) . 


2.0 

3.0  4.0  5.0 

in^  (3.0  Hz) 

Figure  2.4a. 

Variable  frequency  magnitude,  in^Cf),  estimates, 
not  corrected  for  noise,  at  fc  = 0.5  Hz  and 

3.0  Hz  for  events  recorded  at  Bluff,  Alaska 

yl 

(BFAK) . The  arrows  indicate  the  principal 
directions  in  which  the  m^ff)  estimates  would 
move  when  noise  corrections  are  applied.  The 
straight  line  on  this,  and  subsequent  mb(f) 

plots,  marks  the  approximate  boundary  between 
earthquake  and  explosion-like  events. 

1 
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Figure  2.4b.  Variable  frequency  magnitude,  m.  (f),  estimates 
not  corrected  for  noise,  at  0.55  Hz  and  3.0  Hz 
for  events  recorded  at  Bluff,  Alaska  (BFAK) . 


(3.5  Hz) 


Figure  2.4c.  Variable  frequency  magnitude,  m^ff),  estimates 
not  corrected  for  noise,  at  0.5  Hz  and  3.5  Hz 
for  events  recorded  at  Bluff,  Alaska  (BFAK) . 


general  the  noise  correction  required  for  events  in  the  explo- 
sion like  population  primarily  affects  the  low  frequency  mag- 
nitude estimate.  For  events  in  the  earthquake  population,  the 
required  noise  correction  affects  either  both  the  high  and  low 
frequency  estimates  or  the  high  frequency  estimate  preferen- 
tially. 

As  described  above,  we  would  expect  the  application  of 
noise  corrections  to  increase  the  separation  of  the  two  event 
populations.  To  test  this  we  computed  mean  noise  levels  for 
all  the  explosion-like  events  (except  number  1;  no  noise  seg- 
ment available)  and  several  of  the  earthquakes  in  Figures  2.4a 
and  2.b.  The  means- were  based  on  envelope  peaks  preceding  the 
signal  arrivals  as  noted  in  Figures  2.2b  and  2.3b.  The  noise 
estimates  were  then  subtracted  from  the  signal  amplitudes  at 
the  corresponding  frequencies  and  noise  corrected  m^ff)  values 
were  recomputed.  These  revised  ( f ) estimates  are  plotted  in 

Figures  2.5a  and  2.5b  along  with  population  boundary  lines. 
Comparing  Figures  2.4  and  2.5  it  is  immediately  apparent  that 
the  separation  of  the  event  populations  has  been  significantly 
increased  by  the  application  of  noise  corrections.  In  addition 
there  is  a definite  reduction  in  the  scatter  of  the  explosion- 
like population. 

Figures  2.6  through  2.11  give  in^ ( f ) results  for  several 
of  the  stations  previously  reported  on  (LASA,  ANMO  and  KAAO) 
as  well  as  four  new  ones  (ILPA,  KSRS,  CTAO  and  CHTO) . As  of 
the  end  of  this  reporting  period  no  noise  corrections  had  been 
applied  to  the  events  in  Figures  2.6  through  2.11.  As  a 
result  the  separation  of  event  populations  is  not  nearly  as 
large  as  that  in  Figures  2.5a  and  2.5b.  Nevertheless,  it  is 
still  possible  to  draw  some  conclusions  about  many  of  the 
events  with  lower  signal-to-noise  ratios.  Before  summarizing 
the  discrimination  results  there  are  several  points  to  be 
noted  in  particular.  First,  comparing  the  results  in  Figures 
2.4  through  2.11,  we  see  that  the  only  well  recorded  events 


(3.0  Hz) 


Figure 


.5a.  Noise  corrected  m^ff)  estimates  for  events 
recorded  at  BFAK.  Compare  this  figure  with 
Figure  2.4a. 
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Figure  2.6.  Variable  frequency  magnitude,  m^f),  estimates 
not  corrected  for  noise,  at  fc  = 0.45  Hz  and 
2.25  Hz  based  on  short-period  P waves  recorded 
at  ANMO  and  LAO.  These  estimates  are  compared 
with  event  populations  (dashed  and  solid  lines) 
previously  studied  at  LASA  (see  Rodi,  et  al. , 
1978). 
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Figure  2.7.  Variable  frequency  magnitude,  mjj(f),  estimates, 
not  corrected  for  noise,  for  events  recorded  at 
the  Iranian  array  (ILPA) . 
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Figure  2.10. 


Variable  frequency  magnitude,  m^(f) , estimates, 
not  corrected  for  noise,  for  events  recorded  at 
the  SRO  station  in  Charters  Towers,  Australia 
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Figure  2.11.  Variable  frequency  magnitude,  m^Cf),  estimates, 
not  corrected  for  noise,  for  events  recorded  at 
the  SRO  station  in  Chaing  Mai,  Thailand  (CHTO). 


32 


that  cross  the  population  boundaries  (i.e.,  the  straight  lines) 
are  events  that  are  thought  to  be  earthquakes.  For  instance, 
event  63,  which  occurs  landward  of  the  trench  off  northern 
Japan  plots  in  with  the  exp lesion- like  events  at  stations 
BFAK,  ANMO  and  KAAO,  but  plots  in  the  earthquake  population 
at  ILPA.  Event  47,  located  in  the  Kuril  Islands,  would  be 
identified  as  an  earthquake  at  LASA,  ANMO,  BFAK  and  ILPA,  but 
explosion- like  at  CTAO  and  CHTO.  On  the  other  hand,  in  no 
case  does  a well  recorded  explosion- like  event  cross  over  into 
the  earthquake  populations  at  any  of  the  stations  analyzed. 


2 . 4 SUMMARY 

During  this  reporting  period,  we  received  short-period 
seismograms  in  digital  format  for  nine  additional  Eurasian 
events,  bringing  the  total  number  of  events  to  28,  and  con- 
tinued with  an  event  identification  study  based  on  variable 
frequency  magnitude  estimates.  Table  2.3  summarizes  the  pre- 
liminary event  identifications  that  have  been  made  to  date 
based  on  results  from  eight  of  the  stations  contributing  data 
to  this  experiment.  The  most  important  thing  to  note  in  this 
table  is  the  fairly  consistent  dichotomy  of  most  of  the  events. 
Of  the  three  events  (47,  49  and  63)  which  do  cross  population 
boundaries,  the  majority  of  stations  predict  that  two  (47  and 
49)  are  earthquakes  and  the  third  (63)  explosion-like.  The 
locations  of  these  three  events  (Figure  2.1)  along  the  Japan 
and  Kuril  arcs  suggests  that  they  are  most  likely  earthquakes. 
Thus,  it  is  clear  that  a multi-station  discriminant  would 
improve  on  the  results  in  Table  2.3.  At  this  point  in  time, 
however,  with  only  28  events  in  the  data  base,  and  fewer  than 
28  available  for  any  one  station  (Table  2.3),  we  cannot  define 
the  population  statistics  that  are  necessary  for  a multi- 
station analysis.  Assuming  the  data  base  continues  to  increase 
at  its  past  rate  we  should  be  able  to  test  a multi-station 


discriminant  at  KAAO,  KSRS,  ILPA  and  the  Alaskan  stations  within 
the  next  reporting  period. 

Our  plans  for  the  futu:*e  are  the  following: 

1.  Continue  to  test  for  discrimination  at  all 
the  stations  participating  in  the  discrimina- 
tion experiment. 

2.  As  seen  in  Figures  2.4a  and  2.4b,  the  appli- 
cation of  noise  corrections  to  the  ( f ) 
estimates  is  essential  for  optimum  discrimina- 
tion. Thus,  during  the  next  time  period  we 
will  concentrate  on  correcting  the  data  at  all 
the  stations  prior  to  testing  for  discrimina- 
tion. 

3.  Test  for  multi-station  discrimination  using 

a larger  event  population  and  noise  corrected 
m^Cf)  estimates. 
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III.  SURFACE  WAVE  SOLUTION  FOR  THE  ELASTIC 
EQUIVALENT  OF  A COMPLEX  AXISYMMETRIC  SOURCE 

3.1  INTRODUCTION 

Rodi,  et  al , (1978)  presented  a technique  for  linking 

numerical  source  calculations  with  analytical  techniques  for 
elastic  wave  propagation.  The  objective  is  to  analytically 
continue  the  source  displacement  field  into  an  elastic  region 
where  motion  due  to  propagating  waves  can  be  computed  ef- 
ficiently. The  method  employs  the  elastodynamic  integral 
representation  theorem  of  Burridge  and  Knopoff  (1964)  which 
relates  the  total  displacement  field  inside  a given  volume  to 
the  values  of  stress  and  displacement  on  the  surface  of  the 
volume.  For  the  specific  case  of  an  axisymmetric  source, 
displacements  outside  the  source  region  can  be  computed  given 
only  stresses  and  displacements  monitored  on  the  edges  of 
a cylinder  enclosing  the  source  (which  can  be  obtained  from 
finite  difference  or  finite  element  calculations)  and  the 
point  force  Green's  functions  integrated  around  the  cylinder 
radius  (which  can  be  obtained  by  analytical  methods) . 

The  mathematical  development  of  the  body  wave  problem 
using  Cagniard-deHoop  inversion  techniques  is  given  in  Rodi, 
et  al . (1978).  Here  we  deal  with  the  case  of  surface  waves 

excited  in  a multilayered  elastic  or  linearly  anelastic 
medium.  Though  an  incomplete  description  of  the  motion  in  a 
multilayered  half-space,  the  surface  wave  contributions 
often  can  be  used  effectively  to  infer  source  parameters. 

At  regional  and  teleseismic  distances  the  part  of  the  observed 
motion  due  to  surface  waves  can  usually  be  isolated  for  study 
(Rodi,  et  al.  (1978)).  In  certain  cases,  dominant  motion  at 
very  close  distances  may  be  approximated  quite  well  with  a 
surface  wave  description  (Swanger  and  Boore  (1978)). 
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The  development  to  follow  is  quite  general;  it  can  be 
applied  to  many  cases  where  body  wave  methods  cannot  be  ap- 
plied without  great  difficulty.  The  elastic  region  need  not 
be  homogeneous  in  the  vicinity  of  the  source  region.  Plane 
layer  boundaries  may  intersect  the  source  region  where  the 
displacements  and  stresses  are  monitored.  In  theory  the 
method  could  be  applied  to  continuously  varying  material 
properties  with  depth  in  the  elastic  region.  For  this  case 
there  is  a practical  limitation  in  that  no  efficient  techni- 
que for  computation  of  surface  wave  dispersion  parameters 
exists  for  media  with  continuously  varying  moduli  with  depth, 
though  recently,  methods  have  been  proposed  (Wiggins,  1976). 
Also,  with  only  minor  adjustments  in  the  source  related 
terms,  the  method  could  be  extended  to  full  wave  solutions 
applied  in  the  Fourier  frequency  domain  (Apsel,  et  al.  (1977)). 

3.2  THEORETICAL  DEVELOPMENT 

Given  monitored  displacements  and  stresses  in  the 
source  region,  the  surface  wave  motion  can  be  computed  if  the 
Green's  functions  appropriate  to  surface  wave  motion  can  be 
applied  directly.  For  the  body  wave  problem,  it  is  most  ef- 
ficient to  compute  the  displacement  motion  due  to  one  cylinder 
azimuth  and  integrate  over  azimuth  numerically.  In  the  sur- 
face wave  problem  we  can  analytically  integrate  the  Green's 
functions  over  azimuth  before  computing  any  of  the  final 
displacement  field.  To  show  this,  we  begin  with  the  Green's 
function  representation  of  the  displacement  field  due  to  a 
point  force.  Expressions  here  may  differ  slightly  from  those 
of  Rodi,  et  al . (1978)  because  of  algebraic  errors  present 
in  that  work  and  also  because  of  differences  in  the  coordinate 
system  chosen  here. 

We  choose  a coordinate  system  with  positive  z down- 
ward and  receiver  at  a horizontal  distance  r from  the  source. 
Given  either  a vertical  force  or  a horizontal  force  in  a 


direction  <f)  from  the  receiver  azimuth  at  depth  zQ/  the 
Green's  tensor  components  associated  with  vertical  and  radial 
displacements  observed  at  z can  be  written  in  the  Fourier 
frequency  domain  as 


4irpcu 


00 

/ kdk  [a?  ■ (k,a))  + A?.(k,u>)l  (3.1) 

no)  J ' L ID  J 


where 


A?  . = — 

ID 


9 dj  (kr) 

-v  |z-zn|  -k2e  ^ — -tt 

a'  O'  dkr 


^ ■— — 


with 


ekv^e  J.^  (kr) 


ekv^  J1(kr) 


< Jo(kr> 


9 k‘j. (kr) 

Jo(kr)  - — hr- 


-ekv^e-1^  J-^kr) 


■)  6_1 


-ekVg  (kr) 


-k  JQ (kr) 


/ 7 2 \ 1/2 

■(*'•?)  • 


e = sgn  (z-zQ) 


and  only  the  real  part  is  retained  in  exponentials  involving 
the  azimuth  <j>.  The  subscripts  i,j  are  interpreted  as 
follows : 

i=l  implies  horizontal  displacement  in  the  0=0 
direction; 

i=2  vertical  displacement; 

j=l  horizontal  force  in  the  0 direction; 

j=2  vertical  force. 
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Since  we  are  dealing  with  an  axisymmetric  source,  we 
have  not  considered  the  tangential  components  of  displace- 
ments. Also,  because  of  the  source  symmetries,  we  need  only 
consider  the  contributions  due  to  Rayleigh  waves  in  calculat- 
ing surface  wave  motion.  It  would  seem  appropriate  to  use 
only  the  P-SV  contributions  of  the  Green's  tensor,  and  not 
the  complete  representation  of  the  components  of  interest  as 
given  above.  Recent  work  by  Herrmann  (1978)  has  shown  that 
the  complete  P-SV  and  SH  components  of  the  motion  contain 
non-propagating  near-field  terms  which  cancel  in  the  displace- 
ment field  only  when  the  two  systems  of  motion  are  considered 
simultaneously . In  the  above  expressions  we  have  eliminated 
these  terms.  The  contributions  expressed  above  can  still  be 
considered  P-SV  type  motion,  but  they  are  not  the  complete 
P-SV  motion  in  a formal  sense. 

To  apply  the  above  Green's  functions  to  the  equations 
of  Rodi,  et  al . (1978)  directly,  we  need  to  evaluate  them  on 
the  edge  of  the  cylinder  surrounding  the  source  region  and 
then  integrate  over  all  azimuths.  When  we  displace  the  point 
forces  from  the  origin  to  the  cylinder  edge,  we  must  modify 
the  angular  dependences  (Figure  3.1).  Given  a radial  force 
at  cylinder  radius  a and  azimuth  <t> , the  effective  azimuth 
to  the  receiver  becomes  Also,  the  radial  sense  of  motion 

observed  by  the  receiver  will  change  with  <t>.  We  are  interested 
in  horizontal  motion  in  a fixed  radial  direction  (e  direc- 
tion)  , not  the  effective  radial  direction  for  a given  <)>  (e 

R 

direction).  It  is  clear  from  Figure  3.1  that  the  desired 
contribution  for  a given  <J>  is  cos*);  times  the  calculated 
radial  displacement.  These  corrections  will  be  quite  small 
at  far-field  distance,  but,  because  of  the  vector  summation 
properties  of  Bessel  functions,  they  simplify  evaluation  of 
the  azimuthal  integrations. 


L. 
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SOURCE 

CYLINDER 


Figure  3.1.  Source-receiver  geometry  for 
on  the  source  cylinder. 
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R = (r^  + a}  - 2racos<J>) 

To  evaluate  these  integrals  we  note  the  following  summation 
theorems  for  Bessel  functions  (Gradshteyn  and  Ryzhik,  1965) : 

00 

Jv(kR)ei''*  = £ Jm(ka)  J^kr)*1”"* 
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where  Cv  is  a Gegenbauer  polynomial.  Expanding  the  terms  in 
Equation  (3.2)  with  azimuth  dependences,  we  have 
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In  all  but  one  of  the  expansions  above,  the  <p  dependence  is 
now  only  in  the  argument  of  a complex  exponential.  When 
integrating  over  <p,  all  but  one  term  in  each  expansion  will 
vanish.  The  term  involving  the  Gegenbauer  polynomial  re- 
quires special  attention.  Given  the  identify 

_l  , sin  (m+1)  <p 

(cos+>  * sln»  f > 

the  azimuthal  integration  is  then 
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(Gradshteyn  and  Ryzhik,  1965) . We  can  write 

2 IT  _ ..  OO 

/J,  (kR)  . , o v 

d*  kR  e = (ka)-(kr)  L*  m J2m(ka)  J2m(kr) 

0 m=0 

(3.3) 

For  the  surface  wave  solution  we  are  interested  in  only  far- 
field  contributions.  The  integral  in  Equation  (3.3)  is  bounded 
by  2ir/k(r-a)  so  when  kr  >>  1 and  r >>  a,  this  term  is  neglig- 
ible. It  should  be  noted  that  even  in  near-field  cases,  evalua- 
tion of  this  term  would  present  no  problems , since  the 
infinite  series  on  the  right  will  converge  quite  rapidly  for 
small  values  of  the  arguments. 

Compiling  the  azimuthally  integrated  terms  we  are 
left  with  the  azimuthally  averaged  Green's  functions  which, 
neglecting  the  near  field  term,  can  be  written 


(3.4) 
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To  find  the  surface  wave  contributions  due  to  a 
source  described  above,  we  must  first  consider  the  complete 
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multilayered  problem.  The  contribution  due  to  a source  em- 
bedded in  a multilayered  elastic  medium  can  be  expressed  by 
a homogeneous  solution  valid  within  each  layer  plus  a source 
contribution  inserted  as  a particular  solution  into  the  layer 
contribution  containing  the  source.  When  the  boundary  condi- 
tions are  satisfied  at  each  layer  boundary,  it  can  be  shown 
(Haskell,  1953)  that  the  solution  to  the  displacement  field 
evaluated  at  the  surface  can  be  written  in  the  form 


U(u) 


kdk 


F (k , u) , Zq  ) 
G (k,ui) 


J_(kr) 
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where  the  function  G(k,w)  is  independent  of  any  source 
contributions.  The  surface  wave  contribution  arises  from 
singularities  in  the  above  integral  occurring  when  the 
function  G(k,oj)  vanishes. 


If  we  are  interested  only  in  the  surface  wave  contri- 
bution, we  need  to  find  the  k,u>  pairs  where  the  above  integral 
is  singular  and  compute  the  residues  given  an  appropriately 
deformed  contour  in  the  complex  k plane.  Harkrider  (1964) 
showed  that  the  residues  can  be  expressed  in  a rather  con- 
venient form.  For  example,  the  Rayleigh  wave  Green's  tensor 
components  for  a point  force  embedded  at  depth  z^  and  re- 
ceiver at  the  surface  can  be  written 
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where 
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k (gj) 
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VV 

X(kn,z) 
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= horizontal  wave  number  associated  with 
the  n-th  Rayleigh  mode  for  a frequency  ui 

= medium  amplitude  response  function 

= real  normalized  horizontal  displacement 
eigenfunction 

= normalized  vertical  displacement  eigen- 
function 

= Hankel  function  of  the  second  kind  of  order 
m. 


For  a given  set  of  frequencies , the  programs  of  Harkrider 
(1970)  can  be  used  to  compute  all  of  the  above  parameters 
for  a given  Rayleigh  wave  mode.  Equation  (3.5)  is  the  sur- 
face wave  equivalent  of  Equation  (3.1). 

Since  we  know  that  the  source  contribution  enters  into 
the  complete  solution  in  a multilayered  medium  as  particular 
solution  only,  we  can  write  down  the  surface  wave  solution 
for  an  azimuthally  integrated  axisymmetric  source  immediately 
(surface  wave  equivalent  of  Equation  (3.4)). 


X(0)  X ( zQ)  J1(ka)  h|2) (kr) 
-iX(zQ)  JL(ka)  Hq2) (kr) 


-iX (0 ) W(zQ)  JQ (ka)  h{2) (kr) 
-W ( z Q ) JQ  (ka)  H(J2)  (kr) 


(3.6) 

If  we  make  the  far-field  approximation  kr  >>  1,  we  can  employ 
an  asymptotic  representation  of  the  Hankel  functions 


H(2)  (kr) 
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The  resulting  far-field  surface  wave  Green's  tensor  components 
are 
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-X(0)X(z0) JL(ka)  iX (0) W( zQ) JQ (ka) 
Xfz^J^ka)  -iW(zQ) JQ (ka) 

(3.7) 

Note  that  the  horizontal  and  vertical  spectral  displacements 
differ  only  by  a factor  -X(0) , the  surface  ellipticity  for 
Rayleigh  waves.  Without  loss  of  generality,  we  can  continue 
by  examining  only  the  vertical  component  of  displacement. 
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To  obtain  the  analytical  continuation  of  motion  from 
a cylindrical  source  region,  we  will  need  the  derivatives  of 
the  Green's  tensor  with  respect  to  zQ  and  the  source  radial 
variable,  which  we  will  now  call  rQ : 


and  using  identities  in  Har<.:  if  ■ 


where  t(Zq)  and  I(Zg)  are  the  - r.  r •. a ire  verti 

cal  stress  eigenfunctions,  reaper * . 

Finally , if  we  are  given  monitored  stresses  and  dis- 
placements on  a cylinder  of  radius  a and  depth  b enclosing 
the  source  region,  we  can  use  the  expressions  in  Rodi,  et  al . 
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(1978)  to  calculate  the  displacement  field  in  the  elastic 
medium.  Since  we  are  dealing  with  Fourier  transformed 
Green's  functions  and  monitored  quantities  here,  all  time  con- 
volutions become  products  in  the  formulas.  Substituting  into 
Equation  (3.2)  of  Rodi,  et^  al . (1978)  we  have 
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The  horizontal  spectral  displacement  Ui_(r,0,o))  is  simply 
-X(0)  Uz(r,0,co)  . 

3 . 3 IMPLEMENTATION 

The  formulation  above  is  being  programmed  and  tested. 
Evaluation  of  Equation  (3.8)  is  straightforward  and  should 
present  few  difficulties.  Except  for  the  monitored  displace- 
ments and  stresses,  all  information  needed  is  standard  out- 
put of  the  programs  of  Harkrider  (1970) . For  many  problems 
of  interest.  Equation  (3.8)  could  be  simplified  a great 
deal.  For  long  period  surface  waves,  where  wave  lengths 
are  much  larger  than  the  dimensions  of  the  source  region,  the 
surface  wave  eigenfunctions  needed  would  essentially  be 
constants.  The  stress  eigenfunctions  E(Zg)  and  t(Zq)  could 
be  made  zero,  W(Zg)  unity,  and  X(Zg)  the  surface  ellipticity. 
The  source  contribution  from  the  depth  integral  would  be 
simply  an  average  of  the  monitored  values  with  depth.  The 
radial  integral  could  be  simplified  using  the  J^kr^)  ~ 1 
and  J^kr^  = krQ/2  for  (krQ)2  <<:  1*  Under  these  assumptions 
and  including  only  the  lowest  order  terms,  the  vertical  dis- 
placement spectrum  could  be  approximated  by 


The  algorithm  given  here  requires  the  spectra  of  dis- 
placements and  stresses  computed  by  numerical  methods.  Prob- 
lems occasionally  may  arise  because  often  the  periods  of 
surface  wave  motion  of  interest  may  be  much  longer  than  the 
duration  of  the  computed  source  displacements.  In  such  cases 
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spectral  values  at  the  periods  of  interest  will  have  to  bo 
interpolated  from  discrete  values  computed  numerically.  The 
accuracy  of  interpolation  of  the  long  periods  may  be  rather 
sensitive  to  whether  or  not  static  values  are  reached  in  the 
source  simulations.  Such  problems  will  probably  have  to  be 
dealt  with  on  a case-by-case  basis. 

■ 
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IV.  SOURCE  CALCULATIONS 


We  have  initiated  an  extensive  review  of  our  one-dimen- 
sional spherically  symmetric  calculations  of  the  ground  motion 
and  source  functions  due  to  nuclear  explosions  in  salt  (the 
SALMON  event) , in  NTS  granite  (PILEDRIVER)  and  Rainier  Mesa 


tunnel  tuff.  The  primary  purpose  of  this  review  is  to  update 
and  improve  our  nonlinear  constitutive  models  while  insuring 
that  they  are  compatible  with  our  two-dimensional  codes.  With 
the  updated  models,  we  plan  to  investigate  such  two-dimensional 
effects  as  nonhydrostatic  in  situ  stress  conditions,  depth  of 
burial  and  spall  and  slapdown  at  the  free  surface. 

The  most  important  modeling  changes  to  date  are  the  im- 
plementation of  good  equations-of-state  for  the  cavity  gases 
to  replace  the  constant  y ideal  gas  treatment,  improvements 
which  make  the  overburden  pressure  treatment  consistent  with 
our  equations  of  state,  and  the  development  of  a better  effec- 
tive stress  law  to  model  the  influence  of  water  saturation  in 
rocks.  Here  we  present  discussion  of  our  latest  one-dimensional 
calculations  for  PILEDRIVER,  SALMON,  and  Rainier  Mesa  (Area  12) 
saturated  tuff,  and  discuss  the  modeling  improvements  made  for 
each.  Table  4.1  gives  the  material  properties  used  for  the  source 
function  calculations.  A complete  discussion  of  the  basic  con- 
stitutive modeling  may  be  found  in  Bache  et^  al.  (1975)  or  Cherry, 
Rimer,  and  Wray  (1975). 

PILEDRIVER 

It  is  well  known  that  the  measured  velocity  profiles  for 
PILEDRIVER,  a 61  kt  event  in  NTS  fractured  granidiorite  (Perret, 
1968) , cannot  be  matched  using  the  triaxial  failure  envelope 
measured  in  the  laboratory.  A reasonable  approximation  to  the 
laboratory  failure  envelope  is  a parabola  of  the  form 


MATERIAL  PROPERTIES  USED  FOR  SOURCE  FUNCTION  CALCULATIONS 
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Old  relaxation  effective  stress  law,  cavity  treatment,  and  overburden  treatment 
Voids  used  for  effective  stress  law  only. 


L. 


where  Y is  the  maximum  stress  difference,  YQ,  Y , and  Pm  are 
constants  of  the  fit  (given  in  Table  4.1),  and 

1/3 

■ 

where  is  the  third  deviatonc  stress  invariant. 

Calculation  310  of  Table  4.1  and  calculation  130  of  Bache 
et  al.  (1975)  were  successful  attempts  to  reproduce  the  free 
field  velocity  data  of  Perret  (1968)  through  the  concept  of 
an  "effective  stress  law."  In  the  highly  nonlinear  near  field 
regime  this  allowed  us  to  use  a much  lower  failure  envelope 
while  retaining  the  laboratory  failure  surface  outside  the 
range  where  the  maximum  stress  is  less  than  Pc.  The  effective 
stress  concept  may  be  simply  stated.  For  pressure  dependent 
failure  envelopes  the  effective  pressure,  P-P  , where  Pw  is 
the  pore  water  pressure,  should  be  used  to  determine  Y rather 
than  the  mean  stress.  For  rocks  containing  air-filled  poros- 
ity <(),  the  maximum  stress  level  Pc  can  reasonably  be  chosen 
as  the  crush  pressure,  the  pressure  at  which  all  air-filled  por- 
osity is  irreversibly  removed,  since  the  pore  water  cam  hardly 
be  expected  to  carry  much  of  the  load  while  voids  remain  open. 

For  PILEDRIVER,  water  is  assumed  present  in  the  existing  frac- 
tures, rather  than  in  the  negligible  porosity. 

The  simple  effective  stress  model  (we  shall  call  this  the 
relaxation  model)  used  in  calculation  310  amd  in  the  DNA  study 
(Bache  et  al.  (1975))  assumed  that  if  the  meam  stress  in  am 
element  ever  exceeded  P , that  material  would  forever  have  a 

G 

zero  effective  stress,  implying  a much  lower  failure  surface. 

The  stress  deviator  would  be  allowed  to  relax  down  to  the  low 
failure  surface  in  some  relaxation  time,  5.  If  the  meam  stress 
never  exceeded  Pc,  the  laboratory  failure  surface  would  be  used. 

This  relaxation  scheme  was  used  for  the  DNA  study  for  NTS  gran- 
ite, Pahute  Mesa  tuffs  and  rhyolites.  Rainier  Mesa  saturated  tuff, 
and  Yucca  Flat  saturated  and  relatively  dry  tuff  and  was  found  to 


give  a consistent  set  of  equivalent  elastic  sources  which 
could  explain  most  of  the  teleseismic  data  from  NTS  events. 
However,  the  relaxation  model  has  several  drawbacks;  (1)  the 
results  depend  on  an  arbitrary  parameter  5 , the  relaxation 
time  which  scales  with  yield.  (2)  The  discontinuous  change  in 
models  at  the  range  of  maximum  pressure  Pc  results  in  enhanced 
tensile  fracturing  in  some  cases.  (3)  For  materials  with 
negligible  air-filled  voids,  such  as  granite  or  rhyolite,  Pc 
itself  is  chosen  somewhat  arbitrarily.  We  avoided  some  of 
these  difficulties  by  calibrating  <S  and  Pc  to  the  PILEDRIVER 
near  field  data  and  using  the  same  <5  for  all  calculations  of 
the  study  (if  the  material  had  no  voids,  the  same  PQ  was  also 
used  throughout) . 

We  have  now  developed  an  effective  stress  law  which,  for  a 
material  having  negligible  air-filled  voids,  has  no  free  para- 
meters. We  relate  the  change  in  air-filled  porosity  directly 
to  the  effective  stress  and  assume  that  below  some  elastic 
pressure,  Pe,  the  effective  stress  is  the  mean  stress.  As  the 
voids  are  crushed,  the  effective  stress  reduces  to  zero  smoothly 
at  the  crush  pressure.  For  PILEDRIVER,  a material  with  negli- 
gible air-filled  porosity,  we  assumed  that  a small  amount  of 
voids  were  present  and  created  a crush  curve.  However,  this 
crush  curve  was  not  used  in  computing  the  material  pressure  but 
was  used  only  for  the  effective  stress  computations  to  give  a 
smooth  transition  between  regions.  The  results  do  not  depend 
much  on  PQ  but  are  quite  dependent  on  Pc.  We  again  calibrated 
our  model  using  the  PILEDRIVER  near  field  velocity  data  and 
cavity  radius  data. 

For  calculation  410  of  Table  4.1  with  P =2  kbars,  the 

c 

computed  cavity  radius  (based  on  a yield  of  20  tons)  was  far 
too  small.  However,  both  calculations  407  and  411  gave  reason- 
able cavity  radii  and  near  field  velocity  profiles  (the  mea- 
sured cavity  radius  scales  to  3.07  meters  at  20  tons).  Figures 
4.1  and  4.2  show  the  calculated  velocity  profiles  compared  to 
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Figure  4.1.  Particle  velocity  versus  time  at  a range  of  668  feet  for  PILEDRIVER 
(61  kt)  calculation  411. 


Figure  4.2.  Particle  velocity  versus  time  at  a range  of  1543  feet  for 
PILEDRIVER  (61  kt)  calculation  411. 


the  measured  velocity  at  two  recording  stations  (Perret,  1968). 
Peak  velocities  are  in  good  agreement  with  the  data.  However, 
overall  the  waveforms  are  not  in  as  good  agreement  with  the 
observations  as  were  the  computed  waveforms  given  by  Bache, 
et  al . , (1975).  Figure  4.3  shows  the  computed  reduced  velocity 
potential  transforms  for  calculations  410,  407  and  411  together 
with  the  results  for  310.  The  study  shows  that  the  higher  the 
crush  pressure,  the  lower  the  RDP  and  the  more  peaked  the  spec- 
trum. These  relationships  are  reasonable  since  a higher  P 

C 

implies  a higher  strength,  therefore  a smaller  cavity  and 
smaller  The  peaking  of  the  spectrum  is  a function  of  the 

amount  of  tensile  cracking,  since  a higher  strength  implies 
more  cracking,  since  a higher  strength  implies  more  cracking, 
the  results  are  consistent  with  intuition.  Either  407  or  411 
is  acceptable  based  on  near  field  data,  cavity  size,  peak 
velocity,  stress,  etc. 

There  are  two  other  major  modeling  differences  between 
310  and  the  more  recent  calculations.  The  newer  calculations 
have  considerably  better  treatment  of  the  hydrostatic  over- 
burden and  of  the  cavity  gases.  In  the  older  calculations, 
the  scalar  overburden  pressure  was  simply  added  to  the  pres- 
sure obtained  from  the  granite  equation  of  state.  For  the 
newer  calculations , the  ambient  rock  was  compressed  initially 
in  the  code  in  order  to  recover  the  overburden  pressure  directly 
from  the  equation  of  state.  This  procedure  is  more  consistent 
and  will  allow  us  to  imput  simply  a depth  dependent  prestress 
for  the  two-dimensional  calculations. 

Finally,  the  newer  calculations  have  a far  better  equa- 
tion of  state  for  the  cavity.  Until  recently  we  placed  the 
device  energy  in  a rock  sphere  having  a mass  of  70  metric  tons 
per  kiloton  of  device  yield  and  calculated  the  pressure  using 
an  ideal  gas  equation  of  state  with  a constant  y of  1.4.  We 
are  now  using  the  quartz  (S^02)  equation  of  state  to  describe 
the  cavity  gases.  This  equation  of  state,  developed  by  Pyatt 
and  Baker  (1978) , models  the  rock  behavior  from  gas  pressures 
of  many  megabars  down  to  pressures  of  several  bars  including 
the  known  phase  changes.  At  present  we  are  still  placing  the 
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Figure  4.3. 


Equivalent  source  functions  for  PILEDRIVER  calcu- 
lations. The  frequency  axis  is  scaled  to  61  kt, 
the  amplitude  axis  to  0.02  kt.  Calculations  410, 
407  and  411  show  the  effect  of  decreasing  Pc. 

Also  shown  is  the  Mueller/Murphy  source  function 
scaled  to  the  PILEDRIVER  yield  and  depth. 
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device  energy  in  a 70  ton  per  kiloton  rock  sphere.  Tests  are 
now  underway  which  examine  the  importance  of  a smaller  initial 
source. 

Also  shown  in  Figure  4.3  is  the  Mueller/Murphy  source 
function  (Mueller  and  Murphy,  1971)  scaled  to  the  PILEDRIVER 
yield  and  depth.  This  source  function  is  based  on  a fit  to 
the  SHOAL  data  and  gives  a reasonable  approximation  to  the 
PILEDRIVER  data  (Murphy,  1977) , though  it  should  be  kept  in 
mind  that  the  measured  source  functions  for  PILEDRIVER  vary 
over  an  order  of  magnitude.  The  main  difference  between  the 
Mueller/Murphy  source  function  and  those  computed  is  that  the 
later  peak  at  higher  frequencies. 

SALMON 

SALMON  was  a 5.3  kt  nuclear  event  detonated  in  the  Tatum 
salt  dome  in  southern  Mississippi  in  1964.  We  describe  here 
our  most  recent  calculations  of  the  SALMON  event  and  indicate 
the  direction  of  our  present  efforts.  Numerous  unsuccessful 
attempts  have  been  made  to  calculate  the  SALMON  reduced  dis- 
placement potential  (RDP) . The  basic  difficulty  is  that  the 
calculated  RDP  for  a salt  failure  envelope  as  measured  by 
triaxial  loading  in  the  laboratory  (Pratt,  1978;  Heard,  et 
al. , 1975)  is  significantly  smaller  than  the  measured  RDP 
(see  Murphy  1977) . If  a lower  failure  envelope  is  used  to 
raise  the  calculated  RDP,  the  calculated  cavity  radius  becomes 
too  large  when  compared  with  drillback  measurements  in  the 
SALMON  cavi ty . 

Figure  4 . 4 shows  the  calculated  source  spectra  for  the 
material  properties  data  shown  in  Table  4.1.  Calculation  252 
uses  constant  y ideal  gas  treatment  for  the  cavity  and  has 
been  reported  previously  by  Bache,  Cherry  and  Mason  (1976). 

The  rest  of  the  calculations  used  cavity  equation  of  state 
for  salt  derived  using  the  EIONX  equation  of  state  (Pyatt 
(1966)).  EIONX  is  a simple  mathematic  model  which  incorporates 
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Figure  4.4.  Equivalent  source  functions  for  SALMON  calculations 
The  frequency  axis  is  scaled  to  5.3  kt,  the  ampli- 
tude axis  to  0.02  kt. 


I 

many  of  the  important  features  of  the  more  complicated  Saha 
equation  models.  The  new  overburden  pressure  treatment  was 
also  used  for  these  calculations.  No  effective  stress  law  was 
used  for  the  dry  salt  dome. 

Also  shown  in  Figure  4.4  is  the  Mueller/Murphy  source 
function  for  SALMON  (Mueller  and  Murphy,  1971) . Since  this 
source  function  is  based  on  a fit  to  the  SALMON  near  field 
observations,  it  is  shown  as  a convenient  representation  for 
the  observed  data. 

Calculation  408  used  the  lowest  possible  failure  surface 
which  we  could  justify  (based  on  the  uniaxial  measurements  of 
Heard,  et  al.  (1975).  The  calculated  RDP  of  7.8  m3  for  20  ton 
yield  (see  Table  4.1)  was  still  considerably  lower  than  the 
measured  value  (approximately  11  m3  for  20  tons) . Yet  the 
calculated  cavity  radius  was  significantly  higher  than  2.8  m, 
the  measured  value  scaled  to  20  tons  yield.  Calculation  414 
represents  our  best  guess  for  a failure  envelope,  based  on  the 
available  triaxial  data,  in  particular  the  data  from  Terra  Tek, 
Incorporated  (Pratt,  1978)  and  gives  a reasonable  cavity  radius, 
peak  stresses  and  velocities.  However,  the  calculated  RDP  was 
more  than  a factor  of  two  low. 

We  plan  next  to  investigate  whether  a two-dimensional 

. 

in  situ  stress  field  (in  uniaxial  strain)  can  influence  the 
RDP.  At  shot  depth  there  is  a vertical  stress  of  approximately 
175  bars  and  a horizontal  stress  of  55  bars.  Calculation  415 
(made  with  a scalar  overburden  pressure  of  55  bars)  gave  a 
slightly  higher  to  and  cavity  radius  than  calculation  414. 

! 

However,  unlike  the  rest  of  the  calculations,  the  spectra  was 

quite  peaked  due  to  a large  amount  of  tensile  cracking.  Due 

to  this  cracking,  the  elastic  radius  R was  more  than  twice 

e 

as  large  as  for  calculation  414. 

If  in  situ  stress  proves  not  to  be  the  explanation  of 
our  low  calculated  RDP,  then  it  is  likely  that  our  constitutive 
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models  for  salt  are  inadequate.  In  that  event,  we  plan  to 
include  viscoelastic  effects  in  constitutive  model,  possibly 
through  a Maxwell  solid  approach. 
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Area  12  Tuff 

The  saturated  tuffs  of  Rainier  Mesa  are  unique  at  NTS 
in  that  they  have  been  characterized  by  an  enormous  number  of 
laboratory  measurements  of  material  properties  data.  Although 
the  variation  in  material  properties  from  event  to  event  is 
considerable,  we  have  been  able  in  the  last  few  years  to  com- 
pile a series  of  "average"  Area  12  tunnel  tuff  material  prop- 
erties data  which  tend  to  be  valid  for  the  more  recent  nuclear 
events.  A series  of  SKIPPER  calculations  have  recently  been 
completed  for  these  average  tuff  properties  (Table  4.1,  cal- 
culation 409)  which  look  at  the  effects  of  initial  cavity  size 
and  of  the  new  effective  stress  law  on  source  spectra. 

We  were  able  to  investigate  the  effect  of  initial  cavity 
size  through  the  use  of  the  CHEST  24  tuff  chemical  equilibrium 
equation  of  state  which  was  developed  by  Laird  (1976)  , and 
which  accurately  models  the  rock  behavior  in  pressure  regimes 
from  tens  of  megabars  down  to  a few  tenths  of  a bar.  This 
tabular  equation  of  state  has  been  created  especially  for  use 
with  hydrodynamic  codes  to  study  nuclear  explosion  phenomenology. 
It  couples  an  elaborate  chemical  equilibrium  treatment  with 
steam  tables  and  bulk  modulus  data.  The  CHEST  equation  of  state 
was  used  to  describe  the  tuff  both  inside  and  outside  of  the 
cavity. 

Figure  4 . 5 shows  source  spectra  for  the  tuff  calculations 
listed  in  Table  4.1.  Calculation  127  from  the  Bache,  at  al . , 
(1975)  study  used  a smaller  overburden  pressure  causing  the 
spectra  peak  to  be  at  a lower  frequency.  This  calculation, 
using  the  old  "relaxation"  effective  stress  law,  gave  a cavity 
radius  of  6.12  meters  at  the  20  ton  yield.  This  is  far  greater 
than  the  3.72  to  4.41  meter  range  for  the  measured  cavity  radii. 
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Figure  4.5. 

Equivalent  source  functions  for  Area  12  tuff. 

The  frequency  axis  is  scaled  to  10  kt,  the 
amplitude  axis  to  0.02  kt.  All  calculations 
but  127  (from  Bache,  et  al. , 1975)  have  iden- 
tical material  properties.  Calculation  409 
(no  effective  stress  law)  and  412  (with  new 
effective  stress  law)  have  70  ton/kt  initial 

cavities.  Calculations  412/  416  and  413 
respectively  show  the  effect  of  smaller 
initial  cavity  size. 

jc 
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Calculation  412  used  CHEST  24,  the  new  effective  stress  law, 
and  the  new  overburden  pressure  treatment.  It  gave  a cavity 
radius  of  4.30  meters,  well  within  the  data  range.  Calcula- 
tion 409  had  no  effective  stress  law,  but  otherwise  had  the 
same  modeling  and  data  as  412.  The  calculated  cavity  radius 
of  3.875  meters  was  well  within  the  band  of  measured  data. 

A comparison  of  the  three  calculations  indicates  that 
the  old  effective  stress  law  results  in  more  peaked  spectra 
than  is  seen  for  the  newer  effective  stress  law.  For  no 
effective  stress  law  (Calculation  409)  , the  spectra  is  not 
peaked  at  all.  Calculation  412  with  the  new  effective  stress 
law  is  in  better  agreement  with  free  field  data. 

Calculations  412,  416  and  413  have  identical  modeling, 
except  for  the  size  of  the  initial  source.  For  412,  the 
device  yield  was  placed  in  an  initial  cavity  with  radius 
equivalent  to  70  metric  tons/kt  of  yield,  for  416  in  20.7 
tons/kt,  and  for  413  in  8.75  tons/kt.  The  results  for  cal- 
culations 416  and  413  differ  only  slightly.  They  both  show 
a cavity  radius  slightly  larger  than  the  measured  values  and 
gave  approximately  25  percent  greater  cavity  volume,  4^,  peak 
spectra,  and  volume  inside  the  elastic  radius  than  did  calcu- 
lation 412.  A careful  analysis  was  made  of  these  rather  sur- 
prising results.  It  was  noted  that  the  calculated  melt  for 
the  70  ton/kt  initial  cavity  extended  no  further  than  the 
initial  Lagrangian  boundary  of  the  cavity.  For  the  smaller 
initial  cavities,  the  melt  radius  (which  we  call  the  final 
cavity  radius)  extended  considerably  further.  The  analysis 
indicated  that,  for  the  70  ton/kt  initial  cavity  (412) , the 
energy  density  input  into  the  cavity  cells  was  just  sufficient 
to  vaporize  these  cells.  For  the  20.7  ton/kt  cavity  (416), 
far  more  energy  was  input  into  the  smaller  cavity  then  needed 
to  vaporize  the  mass  present.  However,  the  extra  energy  did 
not  vaporize  further  cells.  Since  the  energy  required  to 
vaporize  the  rock  is  wasted  energy  that  could  otherwise  go 
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into  driving  the  shock  wave,  the  smaller  initial  cavity  drives 
a stronger  shock  wave  and  therefore  gives  a larger  final  cavity 
and  RDP.  The  very  small  initial  cavity  (8.75  tons/kt,  calcula- 
tion 413)  has  a sufficiently  high  energy  density  to  vaporize 
rock  out  almost  to  the  same  radius  as  for  calculation  416.  Thus 
the  results  are  quite  similar. 


The  question  that  arises  is  which  calculational  pro- 
cedure is  correct?  Clearly,  hydrodynamic  codes  do  not  take 
into  account  some  basic  physical  processes  such  as  thermal 
conductivity  or  flow  of  vaporized  water  through  the  rock  mass. 
Thus  it  becomes  difficult  to  compute  both  the  correct  vapori- 
zation radius  and  the  melt  radius  using  these  codes  alone 
(equilibrium  procedures  have  been  developed  which  take  the 
late-time  code  output  and  determine  the  true  melt  by  mixing 
the  reserve  energy  of  the  cavity  with  rock  mass  outside) . 
Further  study  is  needed  in  order  to  resolve  this  problem. 
Meanwhile  we  will  continue  to  use  the  70  ton/kt  rock  gas  model 
for  all  calculations.  (For  the  above  calculations,  this  means 
using  412  which  is  in  better  agreement  with  the  measured  cavity 
radius  data) . Any  small  underestimate  of  the  source  spectra 
which  may  be  inherent  in  this  procedure  is  likely  to  occur  in 
source  calculations  in  all  materials. 
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Summary 


We  have  added  the  following  features  to  our  nonlinear 
constitutive  models: 

1.  A new  effective  stress  law  which  scales 
with  device  yield. 

2.  An  improved  treatment  of  overburden  pressure 
which  is  consistent  with  our  equations  of 
state. 

3.  Better  cavity  equations  of  state  for  granite, 
salt  and  tuff. 

Source  calculations  have  been  made  with  these  improved  models 
for  NTS  granite  (PILEDRIVER) , salt  (SALMON)  and  for  average 
NTS  Area  12  tunnel  tuff  material  properties.  Thes  Icula- 
tions  lead  to  the  following  conclusions. 

1.  We  are  unable  to  calculate  the  SALMON  RDP 
in  one-dimension  with  our  present  models 
using  laboratory  data  for  the  failure 
envelope.  The  calculations  give  the  mea- 
sured cavity  radius  but  too  small  an  RDP. 

2.  We  can  calculate  the  PILEDRIVER  RDP  with 
these  models  and  obtain  reasonable  agree- 
ment with  near  field  velocity  data  and 
with  the  measured  cavity  radius. 

3.  Using  the  new  effective  stress  law,  we 
can  calculate  an  RDP  for  average  NTS 
Area  12  tuff  material  properties  without 
any  free  parameters  and  obtain  good  agree- 
ment with  measured  cavity  radius.  We  were 
unable  to  obtain  good  agreement  with  cavity 
radius  data  in  our  earlier  studies. 


Future  work  will  be  directed  toward  the  following  areas: 

1.  SALMON:  We  will  attempt  to  calculate  this 
event  in  two-dimensions  to  study  the  influence 
of  in  situ  stress  on  the  RDP.  We  may  also 
introduce  viscoelasticity  into  our  constitutive 
models. 

2.  PILEDRJVER:  We  do  not  plan  any  further  one- 
dimensional calculations  to  calibrate  our 
models  further  to  obtain  better  agreement 
with  the  Perret  velocity  histories.  Our 
efforts  will  be  directed  toward  two- 
dimensional  calculations  to  examine  the 
effect  of  depth  of  burial  and  spall  on 

the  far-field  signals. 

3.  Cavity  Modeling:  The  calculations  for  tuff 

have  raised  the  question  of  how  best  to  ► 

initialize  our  calculations.  We  will 
examine  this  question  further. 


I 
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